A compressible two-fluid model for the flnite 
volume simulation of violent aerated flows. 
Analytical properties and numerical results 



Frederic Dias 

^LRC MESO, ENS Cachan, CEA DAM, 61 Av. du President Wilson, F-94230 

Cachan, FRANCE 

Denys Dutykh^ 

^LRC MESO, ENS Cachan, CEA DAM, 61 Av. du President Wilson, F-94230 

Cachan, FRANCE 

Jean-Michel Ghidaglia*^ 

""LRC MESO, ENS Cachan, CEA DAM, 61 Av. du President Wilson, F-94230 

Cachan, FRANCE 



Abstract 

In the study of ocean wave impact on structures, one often uses Froude scaling 
since the dominant force is gravity. However the presence of trapped or entrained 
air in the water can significantly modify wave impacts. When air is entrained in 
water in the form of small bubbles, the acoustic properties in the water change 
dramatically and for example the speed of sound in the mixture is much smaller 
than in pure water, and even smaller than in pure air. While some work has been 
done to study small-amplitude disturbances in such mixtures, little work has been 
done on large disturbances in air-water mixtures. We propose a basic two-fluid 
model in which both fluids share the same velocities. It is shown that this model 
can successfully mimic water wave impacts on coastal structures. Even though this 
is a model without interface, waves can occur. Their dispersion relation is discussed 
and the formal limit of pure phases (interfacial waves) is considered. The governing 
equations are discretized by a second-order flnite volume method. Numerical results 
are presented. It is shown that this basic model can be used to study violent aerated 
flows, especially by providing fast qualitative estimates. 
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1 Introduction 



One of the challenges in Computational Fluid Dynamics (CFD) is to deter- 
mine efforts exerted by waves on structures, especially coastal structures. The 
flows associated with wave impact can be quite complicated. In particular, 
wave breaking can lead to flows that cannot be described by models like e.g. 
the free-surface Euler or Navier-Stokes equations. In a free-surface model, the 
boundary between the gas (air) and the liquid (water) is a surface. The liq- 
uid flow is assumed to be incompressible, while the gas is represented by a 
medium, above the liquid, in which the pressure is constant (the atmospheric 
pressure in general). Such a description is known to be valid for calculating 
the propagation in the open sea of waves with moderate amplitude, which do 
not break. Clearly it is not satisfactory when waves either break or hit coastal 
structures like offshore platforms, jetties, piers, breakwaters, etc. 

Our goal here is to investigate a relatively simple two-fluid model that can 
handle breaking waves. It belongs to the family of averaged models, in the 
sense that even though the two fluids under consideration are not miscible, 
there exists a length scale e such that each averaging volume (of size e^) con- 
tains representative samples of each of the fluids. Once the averaging process is 
performed, it is assumed that the two fluids share, locally, the same pressure, 
temperature and velocity. Such models are called homogeneous models in the 
literature. They can be seen as limiting cases of more general two-fluid models 
where the fluids could have different temperatures and velocities [9]. Let us 
explain why it can be assumed here that both fluids share the same temper- 
atures and velocities. There are relaxation mechanisms that indeed tend to 
locally equalize these two quantities. Concerning temperatures, these are dif- 
fusion processes and provided no phenomenon is about to produce very strong 
gradients of temperature between the two fluids like e.g a. nuclear reaction in 
one of the two fluids, one can assume that the time scale on which diffusion 
acts is much smaller than the time scale on which the flow is averaged. Sim- 
ilarly, concerning the velocities, drag forces tend to locally equalize the two 
velocities. Hence for flows in which the mean convection velocity is moder- 
ate (a scale of time is built with the mean convection velocity and a typical 
length scale) we are in the case where this time scale is much smaller than the 
time scale on which velocities are equalized through drag forces. Hence, in the 
present model, the partial differential equations, which express conservation 
of mass (1 per fluid), balance of momentum and total energy, read as follows: 
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where the superscripts ± are used to denote hquid and gas respectively. Hence 
a"*" and a~ denote the volume fraction of liquid and gas, respectively, and 
satisfy the condition + a" = 1. We denote by p"^, u, p, e respectively the 
density of each phase, the velocity, the pressure, the specific internal energy, 
g is the acceleration due to gravity, p := a~^p~^ + a~p~ is the total density, 
E = e + is the specific total energy, H := E + p/p is the specific total 
enthalpy. In order to close the system, we assume that the pressure p is given 
as a function of three parameters, namely a = — a~ , p and e: 



We shall discuss in Section |2] how such a function V is determined once the 
two equations of state p = V^{p^,e^) are known. Equations ([I])-© form a 
closed system that we shall use in order to simulate aerated flows. 

The main purpose of this paper is to promote a general point of view, which 
may be useful for various applications in ocean, offshore, coastal and arctic 
engineering. One can say that the late Howell Peregrine was the first to make 
use of this approach. The influence of the presence of air in wave impacts is a 
difficult topic. While it is usually thought that the presence of air softens the 
impact pressures, recent results show that the cushioning effect due to aeration 
via the increased compressibility of the air-water mixture is not necessarily a 
dominant effect |3]. First of all, air may become trapped or entrained in the 
water in different ways, for example as a single bubble trapped against a wall, 
or as a column or cloud of small bubbles. In addition, it is not clear which 
quantity is the most appropriate to measure impacts. For example some re- 
searchers pay more attention to the pressure impulse than to pressure peaks. 
The pressure impulse is defined as the integral of pressure over the short du- 
ration of impact. A long time ago, Bagnold [1] noticed that the maximum 
pressure and impact duration differed from one identical wave impact to the 
next, even in carefully controlled laboratory experiments, while the pressure 
impulse appears to be more repeatable. For sure, the simple one-fluid models 
which are commonly used for examining the peak impacts are no longer ap- 
propriate in the presence of air. There are few studies dealing with two-fluid 
models. An exception is the work by Peregrine and his collaborators. Wood 
et al. [13] used the pressure impulse approach to model a trapped air pocket. 
Peregrine & Thiais [10] examined the effect of entrained air on a particular 
kind of violent water wave impact by considering a filling flow. Bullock et al. 
[1] found pressure reductions when comparing wave impact between fresh and 
salt water where, due to the different properties of the bubbles in the two 
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fluids, the aeration levels being much higher in salt water than in fresh water. 
Bredmose [2] recently performed numerical experiments on a two-fluid system 
which has similarities with the one we will use below. 

The novelty of the present paper is not the finite volume method used below 
but rather the modelling of two-fluid flows. Since the model described below 
involves neither the tracking nor the capture of a free surface, its integration 
is much less costly from the computational point of view. We have chosen to 
report here on the stiffest case. Should the viscosity effects become important, 
they can be taken into account via e.g. a fractional step method. In fact, when 
viscous effects are important, the flow is easier to capture from the numerical 
point of view. 

The paper is organized as follows. Section [2] provides an analytical study of 
the model. In Section [3] we show that in the limit where the function a is 
identically equal to —1 or 1, equations converge to the equations for 

the classical free-surface flow problems. We also study the dispersion relation 
of the new two-fluid model. Section |4] deals with the numerical discretization of 
this model via a finite volume method. Section [5] is devoted to the presentation 
of the results of various numerical simulations. Finally a conclusion ends the 
paper. 



2 Analytical study of the model 

2. 1 The extended equation of state 

It is shown in this section how to determine the function V{a, p, e) in Eq. (I5l) 
once the two equations of state p = e^) are known. We call Eq. ([5]) an 

extended EOS, since V{—l,p,e) =V~{p,e) and 7^(1, p,e) =V^{p,e), where 

p±=p±(p±,e±), T± = r±(p±,e±), (6) 

are the EOS of each fluid. We will use the following prototypical example in 
this paper. Assume that the fluid denoted by the superscript — is an ideal gas: 

p- = (7- - l)p-e-, = CyT~, (7) 

while the fluid denoted by the superscript + obeys to the stiffened gas law 
(Tait's law) [8]: 

p+ + vr+ = (7+-l)p+e+, e+ = C+T+ + 4^, (8) 
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where 7^, Cy , and tt"*" are constants. For example, pure water is well described 
in the vicinity of the normal conditions by taking 7"*" = 7 and 7r+ = 2.1 x 10^ 
Pa. 

Let us now return to the general case. In order to find the function P, there 
are three given quantities: a G [— 1, 1] , p > and e > . Then one solves for 
the four unknowns , the following system of four nonlinear equations: 

(l + a)p+ + (l-a)p-=2p, (9) 

(l + a)p+e+ + (l-a)p~e~=2pe, (10) 

P+(p+,e+)-P-(p-,e-)=0, (11) 

r+(p+,e+)-T"(p-,e-)=0. (12) 

For given values of the pressure p > and the temperature T > 0, we denote 
by IZ'^^p.T) and £^{p,T) the solutions (p^,e^) to: 

P±(p±,e±)=p, r±(p±,e±)=T, (13) 

and then: 



1 -I- ry 1 — cy 

P = -^n^p,T) + ^7^-(p,T) , (14) 

pe = ^7^+(p, T) £+(p, T) + ^^7^"(p, T) £~ip, T) . (15) 

Finally the inversion of this system of equations leads to p = V{a,p,e) and 
T = T{a,p,e). 



Concerning the prototypical case, the following generalization of ([7]) is consid- 
ered: 

p- + 71- = (7- - l)p-e-, e" = CyT- + . (16) 

7 p 

Introducing 7 and vr defined by 

2 1+a I- a 



7(a) — 1 7+ — 1 7~ — 1 
2 Tiia) 1 + a , 1 — a 



(17) 



7(a) — 1 7+ — 1 7— 1 
Eq. ([H]) and ([IS]) then lead to 

P(a,p,e) = (7(a) - l)pe - 7r(a) , (19) 
X pe-(A+(a)7r+ + A-(a)7r-) 
pCy(a) 
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where 



One can easily check that one recovers the equations of state for each fluid in 
the limits a ^ ±1. 



2.2 An hyperbolic system of conservation laws 

In this section, we assume that the system of equations is solved in M?, having 
in mind the numerical computations performed below. However the extension 
to 3D is trivial. The system can be written as 

^ + V-^(w) = 5(w), (23) 

where 

w = := (a+p+,a~p~, pui, pwa, pE) , (24) 

and, for every G M^, 

jF(w) ■ n = {a^p~^u- n, a" p~u - n, pu-nui +pni, pu-nu2 +pn2, pHu-n) , (25) 

S{w) = {0,0, pgi,pg2,pg- u) . (26) 
The Jacobian matrix A(w) • n is defined by 



In order to compute A(w) ■ n, one writes Eq. ( I25ll for J^(w) ■ n in terms of w 
and p: 

J^{w)-n = [wi ■ ,W2 ■ ,W3 ■ \-pni, 

V Wi + W2 Wi + W2 Wi + W2 



W4- 



W3ni + W4n2 , W3ni+W4n2\ 



+ pn„(^5 + p) ^ ^ : ^ ^ . (28) 
W1+W2 / 



Wi +W2 Wi + W2 

The Jacobian matrix (12 7p then has the following expression: 
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A(w) ■ n 



-Un^^ ^^ni ^^n2 

II. p p p ^ p ^ 



where Un = u- n. 

Let us now compute the five derivatives dp/dwi. A systematic way of doing it 
is to introduce a set of five independent physical variables and here we shall 
take: 

= a, V^2=P, ^3 = T, ip4 = Ui, ip5=U2- (29) 

The expressions of the w-s in terms of the (f'jS are algebraic and explicit. 
Hence the Jacobian matrix dwi/dipj can be easily computed. Since dipj/dwi 
is its inverse matrix, one finds easily with the help of a computer algebra 
program that 



9p , 2 



^1 + ^2) + " P X , (30) 



dwi 

^~\ul + ul) + a^pV, (31) 



dw2 2 



* .(r-i)„,. ^^-ir-i)u,, #^ = r-i, (32) 



dws dwi dw5 

where 

,..2^JfL ±Mr_ ,.,.-^0. (33, 

(cf)^-^^7^(7^-l)T=^^^^, (34) 

r-l^(7(a)-l) (35) 

In Eq. (l35l) . we have introduced the speed of sound of the mixture c^, defined 
by 

_1_ ^ (1 + a)7+ (l-«)7- _ ^ 

pc2 2p+(c+)2 ^ 2p-(c7)^ pa2 ' ^ ^ 

with 

2 _ (l + aK(c+)^ (l-a)p-(c.-)^ 
= 2(-,+ -l) + 2h--l) (3^' 
Then one can show that the Jacobian matrix A(w) ■ n has three distinct 
eigenvalues: 

Xl=Un- Cs, A2,3,4 = Un, A5 = M„ + C^, (38) 
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and is diagonalizable on M. The expression of a set of eigenvectors can be 
obtained by using a computer algebra program. 

Remark 1 //7r+ = and n~ = 0, then = and 



7(o)-l ' 



Remark 2 The left hand side of [3^ is positive since pa^ is bounded from 

Mow by + (^-y-')^ 



2.3 Evolution equations for the physical variables 



The system of conservation laws ([T])-(j4l) can be transformed into a set of 
evolution equations for the physical variables. Let us introduce the entropy 
function s{x,t) defined by (compare with Eq. ( ITOi) ) 



2p s = (1 + a)p+s+ + (1 - a)p's' 



Proposition 1 Continuous solutions to (12^-0) satisfy 

Ut + u- Vu+-Vp = g, (39) 
P 

Pt + u-\/p + pclW ■u = 0, (40) 

at + u- Va + (l-a^)(5V-n = 0, (41) 

St + u-Vs = 0, (42) 

where is given by and 6 is given by 



^ l pcg(7 7r+-7+7r ) 

-2 p^p~{cmcjy ■ ^ > 

Remark 3 For pure fluids (a = ^l), Eq. longer relevant and 6 

is not needed. One can check that the speed of sound Cs is then equal to the 
expected speed of sound (of or c~ ) for pure fluids. 

The balance of entropy (l42l) comes from the balance 

ips)t + V ■ (psu) = 0. (44) 

Adding together Eqs ([1]) and ([2]) leads to 

Pt + V- (pu) = 0. (45) 

Combining Eqs (jH]) and leads to Eq. 
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Remark 4 Subtracting Eq. (QP from ^ leads to 

{px)t + V ■ {pxu) = , with X = • (46) 

P 

In the case of smooth solutions, we obtain that 

+ M ■ Vx = , 

which is an alternative to Eq. 

3 Properties of the model 

3.1 Basic state 

From now on, we denote the set of equations by (E). In order to study 

small perturbations around basic smooth and stationary solutions, it is more 
convenient to use the general set of equations (E) rewritten in the physical 
variables a, u, p and s: see Eq. f l5U]l - fH2]) . 

The two-fluid model (E) describes the evolution of mixtures. It can be used 
for example to study waves along a diffuse interface between a gas and a 
liquid. In order to find the dispersion relation for such waves, one first looks 
for rest states. There is an infinity of such rest states. Then, the governing 
equations are linearized around a special class of rest states. Here the situation 
is considerably complicated by the fact that the stationary solutions are not 
uniform in space. Thus, we come up with a linear system of partial differential 
equations which have non-constant coefficients. 

The steady state will be denoted by a^, p"^, p, u and s. The special class of so- 
lutions we are looking for are motionless, uniform in the horizontal coordinates 
and continuously stratified in the vertical direction: 

cr^ = a^{z), p^ = p^{z), p = p{z), u = 0, s = s{z). 



In this case, Eqs fl39l) - p2l) become 

dp 
dz 



-P{z)9 . (47) 



In the case where one makes the assumption that the mixture density is con- 
stant, 

P(^) = 9(1+ ^)P^ + o (1 - ^)P" = Po, 
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where po is some positive constant, the solution to (H7I) is 

Piz) =Po- pogz, (48) 

where po is a constant. 

Equation pHjl combined with (|T9|) leads to 

Po - ^05-2^ = (7(«) - l)Po e - 7r(a) , (49) 

and there are infinitely many solutions. Since we have imposed p{z) = pq, 
it makes sense to also impose a temperature which does not depend on z: 
T_[z) = Tq. Otherwise, there would be some motion due to convection {u ^ 0). 
Thus, Eq. f l20|) leads to the equation 

e = Cv{a)To + — — , (50) 

Po 

which can be used to obtain an equation for a when combined with Eq. fj^9|) . 



3.2 Linearized two- fluid equations 



In this section we linearize the governing equations around the basic steady 
state derived in the previous section. We write down the following perturbation 
of the stationary solution: 

a = a{z)+2/3+. . . , p = p{z)+q+. . . , u = 0+v+..., s = s{z)+cr+. . . , 

where the vector v has the three components (yi,V2,Vs). 

From Eqs (l40l)-( l42|) . it is straightforward to check that q, (3 and a satisfy the 
equations 



dq 
di 



+ v-Vp{z) + PoqIV ■v = 0, 



da 
'di 



+ v ■Vs{z) = Q . 



(51) 
(52) 
(53) 



In order to obtain the velocity perturbation equation, one needs to compute 
carefully the density perturbation. Using that 

l7=(^l (recall that ^ = 0), 
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P=\{1+ a{z) + 2(3 + .. .)n+ipiz) + g + . . . , To) 

+^(1 - a{z) -2(3 + .. .)n-{p{z) + q + 



.To) 



-il + a{z) + 2p + ...) ip^ 



+ ... 



;(l-a(z)-2/3 + 



7 Q 



Po + (3{p^ - P~) + - 



1 / (1 + a)7+ (1 -a)Y 



Having computed p, it is easy to write the equation for v from Eq. fl39!) : 



dv 1 „ 

dt po 



Po 



2poV (ctr 



(54) 



3.3 Dispersion relation 



From now on, the following notation will be used: 



aiz):=pocl, b{z) := -{1 - a^)6, 

c{z):== — , d{z) ■■= — [ . , .2 + / .2 



(55) 



3.3.1 Case without gravity 

Consider first the simple case where the acceleration due to gravity is absent: 
^ = 0. It represents a major simplification since in this case we recover from 
fl51l) - fl54l) a system of partial differential equations with constant coefficients: 



dq 
di 



+ Pocy-v = 0, 



da 
'dt 



0, 



dv 1 „ 

+ _Vg = 0. 
at Po 
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This system can be written in abstract form as 



where w = (g, cr, v). We look for plane wave solutions: 

w(f, t) = woe*^'''^"'^*^ f = (xi,X2,a;3) with X3 = 2, k = {ki, k2, k-s). 
Substituting this ansatz into equation (ISBj) yields 

3 

B(A^)wo = uwq , with M{k) = ^ hMi 

i=l 

or 



^0 aofci 


ao^2 




60^1 




^0^3 











^00 

PO 








^00 

PO 








y^OO 





J 



bo 



It means that u is an eigenvalue and wq the corresponding eigenvector of the 
matrix B^. 



The dispersion relation of the system (l56l) is given by its characteristic poly- 
nomial 



4 / 2 '^0 1^12 

U! [u \k\ 



or co^ (u^ - cl\k\^ 



0. 



Note that here we have not considered any boundary conditions and that the 
vertical direction does not play any particular role. This is why we have been 
looking for perturbations which are periodic in all three directions. In fact 
there is no dispersion in the acoustic waves we have found. 

Remark 5 The computations performed above are in full agreement with the 
computations ^E^-^^El^)- The matrix B is similar to A(0). Note that we ob- 
tained here 6 eigenvalues as opposed to 5 in Section 2.2. This is only due to 
the fact that we performed the computations in 3D in this section as opposed 
to 2D in Section 2.2. 
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3.3.2 General case with gravity 



Let us now consider the general situation where g is different from zero. We 
drop the equation for a since (jSID^dSSD is a strictly "triangular" linear system 
of PDEs. As above we look for periodic solutions of the form 





















[v] 




v) 



z e 



i{k-x—u)t) 



(57) 



In this case x and k have only horizontal components: 

x = (a;i,X2), k= {ki,k2), 

and V = {vi,V2,w). The main difference with the previous case is that the 
amplitude now depends on the vertical coordinate z. It makes the analysis 
more complicated. 

Substituting the expression flFTI) into Eqs flFIl) - flM|) except for flS51) yields the 
following system of ordinary differential equations: 



-lujq - pogw + a[z) [ik ■ Vi2 + — 
V dz 

■ n da ^ , / \ / -r dw 
-lujjJ + —w + b{z) ik ■ vi^2 + 



dz 



dz 



1 dq 



0, 
0, 



-iuvi 2 + ik — = 0, 
Po 



—iujw H — + c{z)g(3 + d{z)gq = . 

Po dz 



The third equation yields the horizontal divergence of v in Fourier space in 
terms of the pressure perturbation q: 



ik ■ Vi 2 = i 



Another algebraic identity can be obtained if we multiply the first equation 
by b{z), the second one by a{z) and subtract them: 

—iujPa{z) + (^^^a{z) + pQgh{z)^w + iujqh{z) = 0. 
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This relation can be used to eliminate, for example, the volume fraction per- 
turbation (3: 

- hiz) ^ f da b(z) \ w 

a{z) \az a[z)/ too 

Thus, the system governing the behaviour of the perturbations w and q is 
given by the two equations 



Po dz V 



a[z\ 



dw _|_ • / I 1^ ^ 
dz \pouJ a{z) 



+ gc[^z 



f da 
\dz 



+ Po9 



Q - 


Po9 . 
w 

a{z) 


b{z) 




w 


a{z)J_ 


iuj 



-U; = 0,(58) 
0.(59) 



Note that this analysis is consistent with the previous one without gravity. 
Indeed, if one takes g = in (l58l) -( !59l) and assumes a periodic dependence in 
z with wavenumber k^, one recovers the previous dispersion relation. 

In order to find the dispersion relation in the general case, these equations 
must completed by boundary conditions, for example 



'ii'bottom — 0; Wtop 



0, 



if the flow occurs between two solid walls. 



One expects to obtain a dispersion relation as a solvability condition for the 
second-order boundary value problem. Unfortunately, at this stage, we cannot 
go much further, except for the particular case where gz/ (cf )^ <^ 1. Then one 
can perform an asymptotic expansion in this small parameter. As a result the 
coefficients are polynomials in z. Restricting to the leading order terms yields 
a system of two ordinary differential equations with coefficients that are affine 
in z. The solution can be obtained in terms of Airy's wave functions Ai and 
Bi. 



3.4 Pure fluid limit 



We show now that the two-fluid model degenerates into the classical water- 
wave equations in the limit of an interface separating two pure fluids. Consider 
the case where a is either 1 or —1. More precisely let 

a := 1 — 2H{z — ri{x,t)) , x=(xi,X2) (60) 

where Ti is the Heaviside step function. Physically this substitution means 
that we consider two pure fluids separated by an interface. It follows that 

a'^a~ = , 1 — a^ = . 
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Substituting the expression (l60l) into the equation (HT!) gives 



r]t + Uh- VhV = w , 
where Uh = (mi,M2) and V/, = {d^^,d^2). 

This equation simply states that there is no mass flux across the interface. 
Incidentally this is no longer true in the case of shock waves. Integrating 
the conservation of momentum equation ([3]) inside a volume moving with 
the flow and enclosing the interface and using the fact there is no mass flux 
across the interface simply leads to the fact that there is no pressure jump 
across the interface. In other words, the pressure is continuous across the 
interface. Integrating the entropy equation inside the same volume enclosing 
the interface and using the fact there is no mass flux across the interface does 
not lead to any new information. 

One can now write Eqs ([2])-(j4]) in each fluid, either in the conservative form 

(p±), + V-(p^n±)=0, (61) 
{p^u^)t + V-{p^u^0u^) + Vp^=p^g, (62) 
{p^s^)t + V-{p^s^u^) = 0, (63) 

or in the more classical form 



pf + {u^ ■ V)p^ + p±V ■ = , (64) 
u^ + iu^ ■V)u^ + ^ = g, (65) 
sf + ■ Vs^ = . (66) 

In these two systems, the superscripts + and — are used for the heavy fluid 
(below the interface) and the light fluid (above the interface) respectively. 

The system of equations we derived is nothing else than the system of a 
discontinuous two-fluid system with an interface located ai z = ri{x,t). Along 
the interface, one has the kinematic and dynamic boundary conditions 



r]t + u^ ■Vhr] = w^ (67) 
p~=p~^ (68) 



This simple computation shows an interesting property of our model: it au- 
tomatically degenerates into a discontinuous two-fluid system where two pure 
compressible phases are separated by an interface. In Appendix [Bl we derive 
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Fig. 1. Dispersion relation (j7ip for the acoustic mode, 
the dispersion relation for this hmit. We choose the approximate rest state 



T] = 0, p =po, u =0, p =po-pQgz, s 



"0 ) 



and assume that the fluid domain is bounded by two horizontal walls located 
at z = ^a^D {D is the total depth of the domain). Let 9 = Pq / Pq and 
introduce 



\ 



v2 ' 



- 1, s: 



\ 



1 - 



\k\Hc- 



(69) 

There is a distinction between three cases. When c~\k\ < uj < cf\k\, one flnds 
the dispersion relation (IB. 151) . which is reproduced here: 



UJ 



9\k\ 



TjtMT-c^o\k\D) - e5+tanh(5>o+|^|Z}; 



:i - ^)5+Tjtan(Tjao |fc|Z})tanh(5>+|fc|D). (70) 



Neglecting the effects due to gravity, the dispersion relation fjTOjl becomes 

T-tan{Tja,\k\D) = ^5+ tanh(5+a+|^|Z}). (71) 

There is an infinite number of solutions because of the presence of the tan- 
gent term ta.n{T~ aQ \k\D) . Since 6 is small in coastal engineering applica- 
tions, we expect T~ ta.n{T~aQ\k\D) to be small. Then either (T~)'^ is small, 
or T~aQ\k\D ^ mr, G Z. 

After some simple calculations, one obtains for (Tj)^ small ui ~ c7|/c|, and for 
T~aQ\k\D ^ mi 



c:\k\, 



• 1 1- n-Kc. 
with hm uj„ — 



|fc|D->0 



Oq D 
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When uj < min(Cg |/c|, cj"|/c|), one finds the dispersion relation (1B.16|) . which is 
reproduced here: 



9\k\ 



S-t8.nh{S-a^\k\D) + eS;!;t8.nh{S;^a^\k\D)^ 

= (1 - e)5j5+tanh(5-ao |^|/^)tanh(5+a+|^|D). (72) 



In the incompressible limit , the speeds of sound cf go to infinity, ^ 1 and 
one recovers the classical dispersion relation for interfacial waves, namely 



(t8inh{a^\k\D) + etanh{a^\k\D)) = (1 -61) tanh(ao tanh(a^|fc|D). 
g\k\ 



Finally, when u > max(c^|A;|, \k\), one finds the dispersion relation flB.17p . 
which is reproduced here: 

^ (Tjtan(Tjao 1^1^) + ^^.+ tan(r>o+|^|Z})) = 

= -{l-9)T+T-tMT-a^\k\D)tan{T^a^\k\D). (73) 
Neglecting the effects due to gravity, the dispersion relation fl73|) becomes 

Tjtan(rj«o|^|Z}) + ^^r+tan(r>o+|^|Z}) =0. (74) 

There is an infinite number of solutions because of the presence of the tangent 
term. Again, since 6 is small in coastal engineering applications, we expect 
Tj tan(T^ao I A;|-D) to be small. Then we are back to the first case. In this third 
case, if we were not prescribing boundary conditions at the bottom and at the 
top, we could look for perturbations which are periodic in all three directions 
with wavenumber kf in the z— direction. Equation ( IB. 131) then gives 

= {cfn\k\^ + {kfn 

In other words, there is a relationship between the two vertical wavenumbers. 
This relationship is 




It is reminiscent of Snell's law, which describes the relationship between the 
angles of incidence and refraction, when referring to waves passing through a 
boundary between two different isotropic media. 
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4 A finite-volume discretization of the model 



Here we describe the discretization of the model ([I])-(!4]) by a standard cell- 
centered finite volume method. The computational domain C M'^ is trian- 
gulated into a set of control volumes: Q = UxerK. We start by integrating 
equation fl23l) on K: 



d 
di 



[ wdn+ y f J^(w) ■ ukl da= f S{w) dVt , (75) 



where Ukl denotes the unit normal vector on K {~\ L pointing into L and 
M{K) = {LeT : area(K H L) ^ 0} . Then, setting 



^K{t) := ] f w(x,t) dn 

vol A JK 



vo\{K) Jk 

we approximate (ITS]) by 

dwK x-^ area(L n J^) ^ , ^ . ^, . ^ 

E — ;;W^$(wx,Wi;fc)=5(w,0, (76) 

where the numerical fiux 

$(wx, wz,; Ukl) ~ ,1 ^ [ ^(w) ■ ukl da , 

area(L fl K) Jkhl 

is explicitly computed by the FVCF formula of Ghidaglia et al. [6]: 

$(v, w; n) = ^ sgn(A„(/i(v, w)))^^ '-^ , 

(77) 

where the Jacobian matrix A„(/i) is defined in ( !27|) . /x(v,w) is an arbitrary 
mean between v and w and sgn(M) is the matrix whose eigenvectors are those 
of M but whose eigenvalues are the signs of the eigenvalues of M. 

So far we have not discussed the case where a control volume K meets the 
boundary of Q. Here we shall only consider the case where this boundary is a 
wall and from the numerical point of view, we only need to find the normal 
flux ■ fi. Since u{x, t) ■ n = for of G dfl , we have 

{^■n)l:s^Qn = (0,0,^6^,^6^12,0), pb := pl^^g^ , 

and following Ghidaglia and Pascal [7j, we can take pt = P + pUnCg, where the 
right-hand side is evaluated in the control volume K. 

Remark 1 In order to turn ( 76 ) into a numerical algorithm, we must at 
least perform time discretization and give an expression for /i(v,w). Since 
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this matter is standard, we do not give the details here but instead refer to 
Dutykh Let us also notice that formula ( [7^ leads to a first-order scheme 
hut in fact we use a MUSCL technique to achieve better accuracy in space IT^ - 



5 Numerical simulations 

5.1 Basic tests 

In order to check the accuracy of our second-order scheme, we first solve 
numerically the scalar linear advection equation 

_ + a • Vw = 0, a e M^ 
ot 

with smooth initial conditions with compact support in order to reduce the 
influence of boundary conditions. It is obvious that this equation will just 
translate the initial form in the direction a. So, we have an analytical solution 
which can be used to quantify the error of the numerical method. On the other 
hand, to measure the convergence rate, we constructed a sequence of refined 
meshes. 

Fig. [2] shows the error of the numerical method in Lqo norm as a function of the 
mesh characteristic size. The slope of the curves represents an approximation 
to the theoretical convergence rate. On this plot, the curve with circles for 
the data points corresponds to the first order upwind scheme while the other 
two correspond to the MUSCL scheme with least-squares and Green-Gauss 
gradient reconstruction procedures respectively. The slope of the curve with 
circles is equal approximatively to 0.97, which means first-order convergence. 
The other two curves have almost the same slope equal to 1.90, thus indicating 
a second-order convergence rate for the MUSCL scheme. Remark that in our 
implementation of the second-order scheme the least-squares reconstruction 
seems to give slightly more accurate results than the Green-Gauss procedure. 

The next figure represents the measured CPU time in seconds, again as a 
function of the mesh size. Obviously, this kind of data is extremely computer 
dependent but the qualitative behaviour is the same on all systems. On Fig. [3] 
one can see that the "fastest" curve is that of the first-order upwind scheme. 
Then we have two almost superimposed curves referring to the second-order 
gradient reconstruction on variables. Here again one can notice that the least- 
squares method is slightly faster than the Green-Gauss procedure. On this 
figure we represented one more curve (the highest one) which corresponds to 
Green-Gauss gradient reconstruction on fluxes (it seems to be very natural in 
the context of FVCF schemes). The numerical tests show that this method is 
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Convergence rate in L norm 
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Fig. 2. Numerical method error in Lc<j norm. 
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Fig. 3. CPU time for different finite volume schemes. 
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□ 



□ 



Mixture density at t = 1 .2 




Fig. 4. Shock tube of Sod: density plot. 



5.2 Falling water column 



The geometry and initial condition for this test case are shown on Fig. [51 Ini- 
tially the velocity field is taken to be zero. The values of the other parameters 
are given in Table lA.ll The mesh used in this computation contained about 
108000 control volumes (in this case they were triangles). The results of this 
simulation are presented on Figures [614TT1 Fig. [12] shows the maximal pressure 
on the right wall ClS db function of time: 



t\ — > max p(x,y,t). 



We performed another computation for a mixture with a"*" = 0.05, a~ = 0.95. 
The pressure is recorded as well and plotted in Fig. [131 One can see that the 
peak value is higher and the impact is more localized in time. 
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1 



0.9 



a 



0.1 
0.9 



a~ 



0.9 
0.1 



0.05 



0.3 0.65 0.7 1 

Fig. 5. Falling water column test case. Geometry and initial condition. 



Fully CDmprsssible homogeneous Ivo phase soh/et 
Mixture density at t = 0,005 




(a) t = 0.005 



Fully compressible homogeneous phase soh/et 
Mixture density at t = 0.060 




(b) t = 0.06 



Fig. 6. Falling water column test case. Initial condition and the beginning of the 
column collapse. 



5.3 Water drop test case 



The geometry and initial condition for this test case are shown on Fig. [TH 
Initially the velocity field is taken to be zero. The values of the other param- 
eters are given in Table lA.li The mesh used in this computation contained 
about 92000 control volumes (again they were triangles). The results of this 
simulation are presented in Figures [T5142T1 In Fig. [22] we plot the maximal 
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Fully c^mprsssible homogmeous Ivo p^ase soh/er 
MiKture density at t = 0,100 



Fully csmprsssible homogeneous t-jvo phase soh/er 
Mixture density alt = 0.125 





(a) t = 0.1 



(b) t = 0.125 



Fig. 7. Falling water column test case. Splash formation due to the interaction with 
the step. 



Fully comprsssible homogeneous phase soh/et 
Mixture density at t = 0.150 




Fully compressible homogeneous \-*fO phase soWet 
Mixture density at t = 0.175 




(a) t = 0.15 (b) t = 0.175 

Fig. 8. Falling water column test case. Water hits the wall. 



Fully compressible tiomogeneous Ivo phase soh/er 
Mixture density at t = 0.200 




Fully comprsssible homogeneous tw phase soWet 
Mixtufie density at t = 0.225 




(a) t = 0.2 (b) t = 0.225 

Fig. 9. Same as Fig. [8] at later times. 

pressure on the bottom as a function of time: 



1 1 — > max p(x, y, t). 

(x,j/)G[0,l]xO^' 
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Fully comprsssible homogmesus t'vo phase soh/er 
Mixture density at t = 0,300 




Fully compressible homogeneous phase soh/et 
MKture density at t = 0.400 




0.6 o.a 



(a) t = 0.3 (b) t = 0.4 

Fig. 10. Falling water column test case. The splash is climbing the wall. 



Fully comprsssible homogeneous Ivo phase soh/et 
Mixture density at t = 0.500 




Fully comprsssible homogeneous t'-vo phase soh/er 
Mixture density at t = 0.675 
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C,2 O.d 



(a) t = 0.5 (b) t = 0.675 

Fig. 11. Falling water column test case. Turbulent mixing process. 
Maximal pressure on the right wall 




0.3 0.4 
t, time 



Fig. 12. Maximal pressure on the right wall as a function of time. Case of a heavy 
gas. 
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2.2 



Maximal pressure on the right wall 
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Q. 



o 



E 



1.8 - 



1.6 



1.2 - 



1.4 



2 - 
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0.3 



0.4 
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0.6 
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t, time 



Fig. 13. Maximal pressure on the right wall as a function of time. Case of a light 
gas. 

It is clear that the pressure exerted on the bottom reaches 2.5po due to the 
drop impact at t ^ 0.16. 

Remark 2 Beginning with Fig. one can see some asymmetry in the solu- 
tion. It is not expected since the initial condition, computational domain and 
forcing term are fully symmetric with respect to the line x = 0.5. This discrep- 
ancy is explained by the use of unstructured meshes in the computation. The 
arbitrariness of the orientation of the triangles introduces small perturbations 
which are sufficient to break the symmetry at the discrete level. 



6 Conclusions 

In this article we have presented a simple mathematical model for simulating 
water wave impacts. Associated to this model, which avoids the costly capture 
of free surfaces, we have built a numerical solver which is: (i) second-order 
accurate on smooth solutions, (ii) stable even for solutions with very strong 
gradients (and solutions with shocks) and (iii) locally exactly conservative 
with respect to the mass of each fluid, momentum and total energy. This last 
property, (iii), which is certainly the most desirable from the physical point of 
view, is an immediate byproduct of our cell-centered flnite volume method. 

We have shown here the good behavior of this framework on simple test cases 
and we are presently working on quantitative comparisons in the context of 
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Fig. 14. Geometry and initial condition for water drop test case. 



Mixltr& densty at t = O.OOS 



Miwue densty at t = 0-0/5 
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0.2 OA 0.6 0.3 



(a) t = 0.005 (b) t = 0.075 

Fig. 15. Water drop test case. Initial configuration and the beginning of the fall, 
real applications. 



A Some technical results 

The constants Cy can be calculated after simple algebraic manipulations of 
equations (I?!), (IE]) and matching with experimental values at normal condi- 
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Mixitredenstyati =0,103 



MiKltre densty at t = 0-1 25 
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0,2 0-4 0,6 O.B 1 



0.2 0.d 0,6 0.3 



(a) t = 0.1 (b) t = 0.125 

Fig. 16. Water drop test case. Drop approaching the bottom of the container. 
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MiMue densty at t = 0.135 
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MiKlue densty at t = 0.1 53 



0,2 0-d 0.6 o.a 



(a) t = 0.135 



0.2 aa 0.6 0.3 i 



(b) t = 0.15 



Fig. 17. Water drop test case. Drop/bottom compressible interaction. 



Mixicre densty at t = 0,1 /5 



MiKtire densty atl =0.200 



IJLJL 



0.2 O.d 0.6 O.E 



0.2 O.a 0.6 0.3 



(a) t = 0.175 (b) t = 0.2 

Fig. 18. Water drop test case. Vertical jets formation. 



tions: 



Po 



(7 -l)Po^o 
+ 7r+ 



(7+ - l)7+Po+To 
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Miwi.re densty at I = 0.225 




MiKltre densty at t = 0-275 



0.2 0.4 o.e o.a i 



(a) t = 0.225 




0.2 o.a 0. 



(b) t = 0.275 



Fig. 19. Water drop test case. Side jets crossing. 



Mixture densty at t = 0325 



MiWLredenstyati =0,550 
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(a) t = 0.325 



(b) t = 0.35 



Fig. 20. Water drop test case. Side jets flowing down the centerline. 



Miwcre densty at t = 0.403 



MiKttre densty at t = 0.450 
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(b) t = 0.45 



Fig. 21. Water drop test case. Central jet reflection from the bottom. 

For example, for an air/ water mixture under normal conditions we have the 
values given in Table lA.li 

The sound velocities in each phase are given by the following formulas: 



(A.l) 
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Maximal pressure on the bottom 




Fig. 22. Water drop test case. Maximum bottom pressm'e as a function of time. 



parameter value 

Po 10^ Pa 

Pq 10^ kg/m^ 

Pq 1.29 kg/m^ 

To 300 K 



1.4 



7+ 



+ 2.1 X 10^ Pa 



166.72 ^ 
C,7 646.0 
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Table A.l 

Values of the parameters for an air/water mixture under normal conditions. The 
rather high value of the acceleration due to gravity does not correspond to any 
physical situation. Nevertheless, this value was chosen in order to accelerate the 
dynamic processes in the test cases. 

B Dispersion relation in the pure fluid limit 



Let us provide the dispersion relation for waves propagating along the inter- 
face in the limit of two superposed compressible heavy fluids. Consider the 
linearization of the equations around the equilibrium state 
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■nix, t) = 0, 

±/- +\ ± ( Qz 
p {x,z,t) =Pq exp ' 



z, t) = 0, 



p^{x, z, t) = {cffpo exp - {cfypo + Po, 



If we assume that gz/{cfy is small, an approximation to the equilibrium state 
is given by 



ri{x,t)=0, (B.l) 

p^{x,z,t)=p^, (B.2) 

i^(f,2,t)=0, (B.3) 

p^{x,z,t)=po- p^gz, (B.4) 

5^ = 4 (B.5) 

We write the following perturbations of the equilibrium state: 

V = V+C+---, p^ = fy^ + g'^ + - ■ ■ , u^ = 0+v^ + ..., = '[y^+q^ + . . . , 

with in addition s"*^ = + a^. 

The linearized equations for both fluids fl64l) - fl66l) read 



^ + Po^V...t + Po^^ = 0, (B.6) 
' 4v.g=^ = o, (B.7) 



dt p^ 
dw'^ 1 dq^ 
dt po 



(B.8) 



The kinematic and dynamic boundary conditions along the interface are 



dC 

^ = w;±(x,0,t), (B.IO) 
q^{x, 0, t) - p+gC = q~{x, 0, t) - p,gC. (B.ll) 
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We perform a classical perturbation analysis: we look for solutions in the form 



i{k-x—u)t) 



that is periodic perturbations with wave number k and angular frequency u. 
One eventually obtains the following second-order ODE for q^{z): 



2/;± 



0. 



(B.13) 



Having g , one can find v^, g and (. Assume a geometry of total depth D 
bounded above and below by rigid walls located at 2; = GqD and z = —a^D 
respectively. The boundary conditions along the horizontal walls are 



w 



[x, =Fa^-D,t) = 0. 



(B.14) 



Satisfying the kinematic and dynamic boundary conditions along the interface 
provides a solvability condition. In other words, the solution to the linearized 
problem provides the dispersion relation uj{k), which shows the presence of 
two kinds of modes: the gravity modes and the acoustic modes. 



The ODE for q"^ shows that the sign of uj"^ — \k\' 



plays an important 



role. There are three cases: both signs are negative, one is positive and one is 
negative, both signs are positive. 



B.l Case where < u; < cf\k\ 



The assumption is based on the fact that the speed of sound is much higher 
than the speed of sound cj. The solution of the linearized problem is 
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q =y4cos T^\k\{z — aQD] 



and 



w 



\k\ 



T-\k\{z-a^D] 



w = lA — —T^ sin 

^ \T-\k\{z-a^D) 



V =A — 3 COS 

1 



g =A- 



cos 



c: 



T-\k\{z-a^D] 



C = A^T- MT-\k\a^D)e''^'-^^-^'\ 
^ Po 



q+ = B cosh S+\k\{z + a^D] 



^i(k-x—uit) 



-tB-^S+ sinh 



S::\k\{z + a^D) 



^i{k-x—Lut) 



'S;^\k\{z + a^D) 



iT = B — -T cosh 
g+ = 5-4-2 cosh tel^K^ + a^D) 



-^i{k-x—Lut) 



^i{k-x—Lut) 



C = sinh(5+| A:|«o+D)e'(^-^^--*). 



Writing that both expressions for ( are equal and satisfying the dynamic condi- 
tion along the interface yields a system of two homogeneous linear equations 
in A and B. The dispersion relation for the acoustic modes is obtained by 
setting the determinant equal to 0: 



9\k\ 

where 9 := Pq/Pq 



(t- iB.Yi{T-a^\k\D) - 9S+ taiih{S^a^\k\D) 

= {1 - 9)S+T-tan{T-a^\k\D)t^nh{S+a+\k\D), (B.15) 



B.2 Case where CO < mm(cj\k\, c^\k\) 



Assume now that u < Cg\k\. The solution of the linearized problem is un- 
changed for the liquid. For the gas it becomes 
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w 



q = A cosh S^\k\{z — aQD 

S~\k\{z-aoD) 



-iA-^—^Su, sinh 



S^\k\{z-aoD) 



V = A — 3 cosh 

g~ = A— cosh S~\k\{z — aQD] 
Co 

C = -A^S-smh{S-\k\a^D). 



Writing that both expressions for ( are equal and satisfying the dynamic con- 
dition along the interface yields a system of two homogeneous linear equations 
in A and B. The dispersion relation for the gravity modes is obtained by set- 
ting the determinant equal to 0: 



9\k\ 



(^S- tanh(S'Jao 1^1^) + ^3+ tanh(S'+a+ 1 k\D) 



;i - e)5j5+tanh(5jao-|^|D)tanh(5+ao+|fc|D). (B.16) 



B.3 Case where uj > ma.x{c'^\k\,Cg\k\) 



Assume now that u > cf\k\. The solution of the linearized problem for the 
liquid becomes 



g+ = 5cos T+\k\{z + a^D] 



\k\ 



w+ = zE^T+sin T+I^Kz + a+D) 
upQ 



v^ = B- 



k 



cos 



g+ = B 



UPQ 
1 

7^ 



cos 



T+\k\{z + a+D) 
T+|^|(z + a+D 



\k\ 



C = -54^Tj-sin(T+|fc|a+/^)e^('^--*). 



Writing that both expressions for C, are equal and satisfying the dynamic condi- 
tion along the interface yields a system of two homogeneous linear equations 
in A and B. The dispersion relation for the acoustic modes is obtained by 
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setting the determinant equal to 0: 

(TJ tan(Tj«o" I ^1 ^) + ^^c.^ tan(T>o+ 1 ^1 1); 



9\k\ 



-(1 - e)T+Tj tan(Tjao|^|D)tan(T>o+|^|D). (B.17) 



C Isentropic flows 



Let us consider the isentropic version of the system of equations ([I])-(j4j). It 
reads 



(aV)t + V- (a+p+M) = 0, (C.l) 
+ V- (a~p~n) = 0, (C.2) 
{pu)t + V ■ {pu®u + 'pl) = pg, (C.3) 

with the equation of state 

p = Vi{a,p), (C.4) 
where the subscript / stands for isentropic. 

One can determine Vj as follows. First consider the two equations 

(l + a)p+ + (l-a)p- = 2p, (C.5) 
Vtip^)-V7ip~) = 0. (C.6) 

Given p > 0, we denote by TZf{p) the solutions p^ to 

Vf{p^)=p. (C.7) 

It follows that 

1 -I- rv 1 — ry 

p = ^ntip) + ^nj{p). (C.8) 

The inversion of this equation leads to p = Vi{a, p), which is the equation of 
state that appears in flC.4p . 

In order to see if one can go further analytically, let us consider the particular 
case of stiffened gases. The equations of state are 

p± + vr± = (7± - l)p±e± , (C.9) 

and 



V ~ + + ■ 
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Proposition 2 Possible entropies for stiffened gases are given by 

(P 



This expression can be easily obtained by integrating the well-known differ- 
ential relation 

Tds = de -|- p d - . 

VP/ 

Remark 6 Other possible entropies for stiffened gases, which differ from 
the previous ones by a constant, are given by 

..^c,tioj -"-:w v ,c„) 



Thus, saying that s is constant boils down to saying that 



where A"^ is a constant. Substituting (1C.12P into the EOS (]C.9I) yields 



± ' ^ = A^ip^r^ , (C.13) 



± 



that is 



Vf{p^) = -^ + A^{p^r\ (C.14) 



Thus 



and the equation which gives Vi{a,p) is 



Even in the special case of two perfect gases where vr^ = 0, this equation is 
in general transcendental. This is to be contrasted with the general case (the 
case with a variable entropy), where V{a, p, e) can be calculated explicitly by 
algebraic equations. 



From now on, we denote the set of equations (]C.l|) - (]C.4p by (E/). In order to 



study small perturbations around basic smooth and stationary solutions, it is 
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more convenient to use the set of isentropic equations (E/) rewritten in the 
physical variables a, u, p. 

These equations are given in the following proposition. 
Proposition 3 The equivalent system to (Ej) in variables a, u, p is 



Ut + u- Vu+-Vp = g, (C.17) 
P 

Pt + u-Vp + pc'j^V -11 = 0, (C.18) 

at + u-Va + {l-a'^)SiV -11=0, (C.19) 

where cj^ and 6j are given by 



4—^^^, r. = l±i^r- + i^r-. (C.20) 

To p+p p 2 2 

r* = ^^^^=7* + -. p = a-/+aV. (C,21) 
p ap^ p 

r+ - 

5i = . (C.22) 

2ro ^ ^ 

Remark 7 In the one-fluid case (take for example = l,a~ = 0,a = 1), 
one finds 

Cs = cis = c+ , (C.23) 
while, in the two- fluid case, Cg 7^ cis- 

The analysis for the dispersion relation is quite similar to the general case. 
The steady state is denoted by a^, p"^, p and u. We look for a special class 
of solutions which are motionless, uniform in the horizontal coordinates and 
continuously stratified in the vertical direction: 

a"^ = cr^iz), p"^ = p^{z), P = p{z), u = 0. 

Again we take p to be constant. One must solve 

{l+a{z))nt{p{z)) + {l-a{z))nj{p{z)) = 2po , (C.24) 

in order to find a{z) and p^. It is easy to see that 

, , '2po-KHp{z))-KJ{p{zj) 



with p{z) given by ( HHl) . 

The analysis is the same as before, except that and 6 are replaced by the 
values for the isentropic case. 
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